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Abstract. The quantum phase transition from the Mott insulator state to the superfluid in 
the Bose-Hubbard model is investigated. We research one, two and three dimensional lattices 
in the truncated Wigner approximation. We compute both kinetic and potential energy and 
they turn out to have a power law behaviour as a function of the transition rate, with the 
power equal to 1/3. The same applies to the total energy in a system with a harmonic trap, 
which is usually present in the experimental set-up. These observations are in agreement with 
the experiment of [8], where such scalings were also observed and the power of the decay was 
numerically close to 1/3. The results confirm the Kibble-Zurek (adiabatic-impulse-adiabatic 
approximation) scenario for this transition. 



1. Introduction 

Excitations resulting from crossing the gapless quantum critical points are a serious problem for 
quantum simulations with ultracold atomic gases or ion traps. There, one would like to prepare 
a simple ground state of a simple initial Hamiltonian and then drive the system adiabatically 
to an interesting final ground state. This general observation has been recently supported by a 
more quantitative theory pQ E] , that is by now confirmed by several numerical studies [3] . The 
theory is a quantum generalization of the classical Kibble-Zurek mechanism (KZM) [HE]. The 
theory predicts that density of excitations or excitation energy decay with a power of quench 
timescale tq. Experiments dedicated to the quantum theory were made in Refs. [U [8]. In 
Ref. [8] ultracold spinless bosonic atoms in a three-dimensional (3D) optical lattice were driven 
from the Mott insulator phase to the superfluid phase across a quantum phase transition. The 
transition was non-adiabatic and the excitation energy was reported to scale algebraically with 
the transition time, and the power was numerically equal to one third. Similarily, adiabaticity 



of loading atoms into an optical lattice [9] or releasing them from the lattice confinement jlO] 
has been questioned recently. The 1/3-scaling reported in the experiment |8j coincides with the 
earlier prediction in Ref. [11] for ID. In [21] it was shown, that the same scalings also hold for 
3D. 

2. The Bose-Hubbard model 

The Bose-Hubbard (BH) model is a lattice model that describes bosonic atoms in a field of an 
optical lattice [U [TUl 121 [13] . Its Hamiltonian assumes the following form 

#BH = ~ J Y1 ( a i fl j + fl j ai ) + 2~ a i a i a i°i + H ^iOiOi , (1) 
<ij> 1 * 

where a\ is the annihilation operator at site i and ai is the corresponding creation operator. Here 
we do not specify the dimensionality of a lattice, which is three dimensional in real systems, but 
can be effectively made two- or one dimensional when some dimensions of the system are very 
small compared to the other. 

The BH Hamiltonian has three terms: (1) the kinetic term, which corresponds to the hopping 
of atoms between the neighbouring sites of the lattice, (2) the on-site interaction term, which is 
responsible for on-site repulsion between the atoms and possibly (3) the interaction with external 
potential, which could be a harmonic trap considered later in the article. Here, the coefficient 
^ in front of the on-site interaction term is a result of our choice of units. 

When the kinetic term dominates over the on-site repulsion, the system is in the superfluid 
regime which is characterized by strong spatial correlations. On the other hand, when the on-site 
interaction term is much larger than the kinetic part, the system is in the Mott Insulator state, 
where there is a definite and constant number of atoms n at each site, 

|MI) = \n,n,n,---) . (2) 

Clearly, between these two phases there is a boundary, crossing which causes the system to 
undergo a phase transition. The phase diagram of (pQ) is depicted in Fig. [H where Mott Insulator 
phase forms the famous lobes, each characterized by a distinct value of atom density n. Those 
lobes are surrounded by a sea of the superfluid phase. 

The relative strength of the two terms is controlled by the coefficient J. To cross the phase 
boundary we vary this quantity in time. In this article, we assume that initially the system is 
in the Mott state and we drive it into the superfluid regime by increasing J(t). This driven 
transition will be called here a quench, and we choose it to be linear in time: 

J=-, (3) 

where tq is the inverse of the quench rate. 

3. The truncated Wigner method 

To investigate the model (pQ), we adopt an approach known as the Truncated Wigner Method, 
where the full quantum dynamics of the original Hamiltonian is approximated by an evolution 
of a statistical ensemble of complex lattice fields fa. This approximation is valid in the regime 
of large atom densities, n> 1, which we now assume. In the Truncated Wigner Method, we 
make an identification «j « y/nfa, a\ « y/nfal- With this identification, the equation of motion 
for ai becomes the nonlinear equation for fa 



id t fa = -jv 2 fa + (\fa\ 2 -i) fa , 



(4) 
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Figure 1. The schematic phase diagram of the system defined by (pQ). The lobes of the Mott 
insulator for J< 1 characterize with a constant number of atoms per site. The surrounding 
superfluid phase has a vanishing energy gap and non-zero correlations in phase between sites. 



which is known as the Gross-Pitaevskii equation (GPE). The quantum evolution of the original 
system (pQ) is represented by the whole set of trajectories of fa and the quantum expectation 
values become statistical averages of relevant observables over the ensemble. The initial condition 
for the evolution becomes the distribution for fa. The initial state ([2j) translates to 

fa(0) = e^, (5) 

where #'s are uniformly distributed over (— tt, tt]. This reflects the uncertainty between the atom 
number and phase: n is well defined and 9 is completely random. 
The kinetic and potential energies take the form: 



E PO t = ^£(W 2 -i) 2 

i 

The behaviour of energy with the rate of the transition will be at the center of our interest. 
4. The impulse-adiabatic scenario 

The next ingredient that we add to our analysis is the impulse-adiabatic scenario [5] . It simplifies 
the picture of the evolution across the critical point. In the vicinity of that point we can treat 
the evolution as impulse, which means that the transition rate is much faster than the pace of 
the change of the system. Far from the critical point, on the contrary, the relaxation rate of 
the system exceeds the driving pace and the evolution is adiabatic. The idea of the impulse- 
adiabatic scenario is to say, that there is a certain value of the parameter, in our case J = J 
where the two rates become comparable, and which separates the impulse stage of the evolution 
from the adiabatic stage. In other words, below J the system undergoes the impulse evolution, 
and above - adiabatic. This idea is summarized in the Fig. [2j 

In our units, the critical value of J is J cr ~ n -2 and tends to zero for large densities, J ~ 0. 



(6) 
(7) 
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Figure 2. The schematic depiction of the impulse-adiabatic approximation to the dynamics of 
the considered model. As long as the relative transition rate, J/J, is larger than the relaxation 
rate, the dynamics is said to be impulse, i.e. the state of the system is frozen in spite of the 
changing Hamiltonian. When the relaxation rate begins to dominate, the state has more and 
more time to follow the Hamiltonian, and we assume that the evolution is adiabatic. 



5. The Josephson regime 

First, we restrict ourselves to the regime, where J <C 1. In this regime the density fluctuations 
are small, \(pi\ 2 ~ 1; the field amplitude does not deviate much from the Mott state's amplitude. 
Hence, it is convenient to parametrize 

A = (1 + /0 , (8) 

with real amplitude correction /; and phase 9\. This gives us equations of motion for f\ and 0{, 
where the initial conditions are f\ = and random 6{. After elimination of /j <C 1 in Eq. (|4|) 
we get the Josephson equations 

|^i = 2J sin^-ft) (9) 

j,n.n. i 

with random initial #i(0) and ^(0) = 0. The summation in Eq. ([9]) extends over the nearest 
neighbours of i. 

The energies defined by Eqs. (JHl EJ) in the approximation ([8]) yield 

E kin « J^VflpV^ , 

i 

^ P ot « 2^^. (10) 

i 

6. Thermalization 

Since the parameter J could be eliminated from Eq. Q by introducing a rescaled time variable 
u = J l l 2 t, the equations have a characteristic time-scale 



~ J- 1 ' 2 



(11) 



This timescale is a relaxation time towards thermal equilibrium in the Josephson regime. In 
one dimension the system has finite correlation length and quickly reaches its equilibrium. The 
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Figure 3. Thermalization in ID after a sudden quench, where J abruptly jumps to its final 
value at t = 0. The short-range correlations quickly reach their equilibrium, whose numerical 
value was predicted in detail in |21j . In each panel, the plots for various J's against time J l l 2 t 
collapse to a single curve, proving, that the thermalization time scales like r ~ J -1 / 2 



different. In 2D the correlations of the system display quasi- long-range order and in 3D they 
are of a long range order. Therefore, the two- or three dimensional system cannot equilibrate 
in finite time. However, short range correlation functions can reach equilibrium in finite time, 
which is shown in Fig 01 The correlations will be discussed in more detail in section PTU1 

Since the kinetic and potential energies are quadratic in their degrees of freedom, Eq. (fTUI) . 
we can apply the equipartition principle and write 

£kin = \TV , 

£ P ot = \TV , (12) 

where we introduced temperature T which is a function of J in the adiabatic quench (tq S> 1) 
and V = L D is the volume of the system, i.e. the number of sites. Eqs. (|12p are valid for the 
Josephson regime, J <C 1. The total energy is, obviously, E = E^m + E po t = TV 

For adiabatic processes, where the system remains in equilibrium for every instant of time, 
both T(t) and J(t) are functions of time t. Hence we can write 



(13) 




Figure 4. Correlations in 3D. The correlation function Cr = c/)*(j) s +ii for several values of the 
rescaled time u = J l l 2 t. The initial state was a result of an evolution from J = to J = 0.01 
with tq = 512. Than it was allowed to thermalize. The expected long range order is not reached 
in finite time. However, the short range correlations do thermalize quickly, and so do the local 
observables, such as the energy. Thus, we can apply the equipartition principle. 



where we assume the energy to be a function of a time dependent temperature and 

± E= <E E , mJ - l = ±LIl. (14) 

dt dt kin dt 2J v ; 

Equating (|13p with (|14p and solving for T(J) we get the equation characterizing the 
thermodynamic process induced by driving J: 

T = A^J, (15) 

which we will call the adiabate equation. Again, it is valid for an adiabatic processes in the 
Josephson regime. 

7. Excitation energy in the Josephson regime 

Having the adiabate equation (fl5j) , we may now attempt to predict the scalings of energies (fTtJj) 
using the impulse-adiabatic approximation introduced in section [5J The transition rate is 

H 

and the relaxation rate of the system 

The cross-over takes place at the instant of time t (corresponding to the parameter value J) 
where these two rates become equal. Thus, 

i- 3/2 - -o 1/2 (18) 



or 

J ~ r Q 2 / 3 . (19) 

The idea behind the impulse-adiabatic scenario is that the initial state remains frozen throughout 
the impulse stage of the evolution. Thus, the adiabatic process begins at J with random phases 
0i, 

~ 2tt 2 ~ 

E kin \j « ^E Vfl i • ™i| e _ random = - 3- JVD . (20) 

i 

The above equation gives the initial temperature for the adiabatic evolution 

4tt 2 ~ 

Tl^—JD. (21) 



Comparing this formula with the adiabate equation allows us to determine the coefficient 

tori! 

3 



1 ^' JD. The adiabate equation becomes 



T ~ —DsfJj ~ v7t q 1/3 , (22) 
3 

where in the last step we used the scaling for J (Eq. (|19p ). As a consequence, we receive the 
algebraic scaling for kinetic and potential energy 

E km ~ JVVq^V , £ pot ~ jVVq^V (23) 

in the Josephson regime. 

The predicted scaling agrees perfectly with the numerical simulation of the full Gross- 
Pitaevskii equation, (j3J), where we also get algebraic scalings with the exponent numerically 
close to 1/3. The calculated E^in and E pot are presented in Fig. [5] and the fitted exponents are 
collected in Table [Q 
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Table 1. The best fits to a in (-E^km/pot 

) ~ r Q a are consistent with a = 1/3. Lattice sizes as 

in Fig. [5l 



8. Excitation energy in the Rabi regime 

The data presented in Fig. [5] suggest, that the Tq scaling extends also to the Rabi regime, 
J> 1. This can also be explained in terms of the impulse-adiabatic scenario. Since, for J> 1, 
the kinetic energy dominates, (E) (E kin ), and (p3|) becomes 



d_ 
~dt 



{E\m 



dJ (gkin) 



(24) 



This has a simple solution (.Ekm) = -SJ", where B is a proportionality constant. The crossover 
from the Josephson to the Rabi regime is at J ~ 1 and thus we can determine the constant B: 
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Figure 5. The numerical evolution of the Eq. leads to the kinetic and potential energies 
decay algebraically with tq. The exponent is numerically equal to — |, see Tab.[TJ in agreement 
with the prediction (|23[) . The lattice sizes were L = 4096, 256, 128 for 1, 2 and 3 dimensions 
respectively. 



— 1/3 

Hence, in the Rabi regime, (E) rj (-Ekm) ~ Jtq V- This is true, as long as the crossover 

2/3 

takes place in the Josephson regime, which according to J ~ Tq , happens for quenches slow 
enough. 

For fast quenches, the impulse-adiabatic crossover takes place in the Rabi regime. The 
thermalization time is set by the strength of the nonlinearity in (JH) (t ~ 1) and, consequently, 



ldJ_ 
J~dt 



T- 1 ~ 1 . 



(26) 



Therefore, t ~ 1. Since at J phases remain random, (-Ekii 



JV. Finally, 



DJV . 



(27) 



For fast quenches, energy does not depend on tq, consistently with the data from Fig. [5j 



9. Excitation energy for the system in a harmonic trap 

In order to make our analysis even more realistic and corresponding to the experimental 
conditions, we place our system in an external harmonic trap, i.e. quadratic potential. This 



amounts to setting the V\ in the Bose-Hubbard Hamiltonian ([I]) to 



Vi = \u?i 2 , (28) 

where to is the frequency of the trap, which sets its width, and we assume i = at the center of 
the lattice. Now, the Gross-Pitaevskii equation becomes 



idtfa = -JV 2 ^ + (|<^| 2 - l) X -uj 2 \ 2 ^ , (29) 



where the initial conditions are again random phases 0\ and the density is distributed according 
to 

l0i(O)| 2 = y(^F-i 2 ), (30) 

for i 2 < B^ F and |</>i(0)| 2 = for i 2 > E^ F . This is an equilibrium solution to ([29]) for J = 0. 
In the above distribution i?TF is a radius (Thomas-Fermi radius) of a sphere inside which the 
distribution is non-zero. We choose uj 2 = p 2 — , so that in the center of the trap initial field 

-TI-tf 

density is equal to 1, |(/>(0)| 2 = 1. This makes comparison with the uniform case easier. 

In order to examine scaling of excitation energy in case of the presence of the harmonic trap, 
we have to subtract the ground state energy from the total energy calculated in the numerical 
simulation. The ground state energy can be computed as a solution to 

which minimizes the energy functional 

= E jjV#V0i + + yi^i j , (32) 

subject to the constraint 

^ \ 4>i\ 2 = const. , (33) 



i.e. the norm should be kept constant throughout the evolution. The minimalization of (|32|) can 
be done by evolving the following equation, 

dt s<k ' 1 ' 

until 0i reaches the steady state, which corresponds to the solution of (]3ip and therefore is the 
ground state. 

The evolving cloud in the trap is initially centered in the form of a Thomas- Fermi distribution, 
but subsequently expands and asymptotically tends to the Gaussian distribution for large J. 
This is depicted in the Fig [3 

The numerical simulation of (I29p indicates, that as in the uniform case, the excitation energy 
decays like Tq with the quench rate. This data is shown in the Fig. [HJ The scaling obtained 
for a uniform system survives also in case of a trap, because the kinetic energy accumulated 
in random phase is much larger than the energy related to the localization of the cloud in the 
center of the trap. 
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Figure 6. Relaxation to a ground state of the energy functional E[(f>i] = 
J2i |«/V0?V0i + ^4>i <j>\ 4>i4>i + ^i 2 0?0i| according to the evolution of the equation ^ = 

— -jjjr-- Three panels represent subsequent field configuration as it converges to the energy 
minimum and the bottom right panel shows the corresponding convergence of the ground state 
energy. The initial configuration was an educated guess: we assumed it is proportional to the 
ground state of the equation in the J> 1, ie. the harmonic oscillator equation. The larger J 
the final ground state has more resemblance to the Gauss function as the harmonic oscillator 
approximation becomes more accurate. 



10. Correlations and Vortices 

In one dimension, the correlation length is finite, i.e. the correlation function is exponential, see 
derivation in [21j, 

C R = -(a\a i ) =e~f , (35) 



with £ depending on J and tq through 



n 



e = % - j 1 ^ 3 ■ (36) 

— 1/3 

With growing tq the correlation length decays like Tq , in agreement with the numerical data, 
presented in Fig.[9l This is not the case, however, in 2 and 3D, where there is a (quasi-)long-range 
order. Nevertheless, due to the finite rate, with which the system is driven, the correlations also 
spread with a finite rate and their range is limited, see Fig. [TUl 

The finite rate of the correlation build-up results in the domain formation and appearance 
of topological vortices. Such vortices, in case of a 2D system, are presented in Fig. [TTJ For fast 
quenches (left panel), the there is little time for the correlations to spread, so the resulting pattern 
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Figure 7. The evolution of the cloud of atoms in a harmonic trap. The figure shows a cross- 
section through the center of a 3D lattice. The initial Thomas-Fermi profile at J = expands 
and for large J (here J = 10) tends to a Gaussian. The simulation for tq = 102.4. 
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Figure 8. The excitation energy (E) — Eqs for the atoms in a trap. The energy scales 
algebraically with tq and the powers are 0.32 for J = 1.0 and 0.31 for J = 0.1 and J = 10.0. 
The initial radius of the Thomas-Fermi profile was R = 10. 



contains small domains and many vortices. When the system is allowed to cool down during 
a slow quench (right panel), the domains are larger and there are fewer, but more pronounced 
vortices. 
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Figure 9. The correlation functions in ID at J = 0.1 are exponential, with the correlation 
length numerically scaling as £ ~ Tq 1 ^ , consistently with (136[) . 
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Figure 10. Correlations at J = 0.1 in 2 and 3D respectively. The correlations have a finite 
range due to a limited pace of their growth. The local degrees of freedom have equilibrated, but 
at large scales the (quasi-)long-range order cannot be reached in finite time. 



11. Conclusion 

The physics beyond the linear quench from the Mott insulator phase to the superfluid turns out 
to be well approximated by the impulse-adiabatic scenario, where the phases remain frozen at 
random values, as for the initial state, but at some point the evolution becomes adiabatic and the 
system is allowed to thermalize. In 2 and 3 dimensions, however, the finite rate of the transition 
does prevent the (quasi-)long-range order to fully develop, and the resulting correlations have 
finite range. This is spectacularly manifested by the formation of vortices in phase. The simple 
picture of the impulse-adiabatic approximation allows for derivation of the algebraic decay of 

— 1/3 

energies with the quench time tq, i.e. E kin / pot ~ Tq . This result is reproduced also when the 
system is subject to the external harmonic potential. This is in agreement with the experimental 
data obtained in [8]. 
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Figure 11. The phase 6 S for a 2D 256 x 256 system. The top panel shows the phase at J = 0.1 
after a relatively fast quench with tq = 1638.4. In this case the domains that formed are small, 

1/3 

as is the resulting correlation length £ ~ Tq . The bottom panel shows the situation after a 
quench with tq = 52428.8. The correlation length is much larger, the domains more prominent 
and the resulting vortices clearly visible (pointed with arrows). 
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